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A key phenomenon related to the Josephson effect is oscillations of different properties of super- 
conducting tunneling junctions with magnetic field. We consider magnetic oscillations of the critical 
current in stacks of intrinsic Josephson junctions, which are realized in mesas fabricated from lay- 
ered high-temperature superconductors. The oscillation behavior is very different from the case of 
a single junction. Depending on the stack lateral size, oscillations may have either the period of 
half flux quantum per junction (wide-stack regime) or one flux quantum per junction (narrow-stack 
regime). We study in detail the crossover between these two regimes. Typical size separating the 
regimes is proportional to magnetic field meaning that the crossover can be driven by the magnetic 
field. In the narrow-stack regime the lattice structure experiences periodic series of phase transitions 
between aligned rectangular configuration and triangular configuration. Triangular configurations 
in this regime are realized only in narrow regions near magnetic-field values corresponding to integer 
number of flux quanta per junction. 



I. INTRODUCTION 

Layered high-temperature superconducting materials, 
such as Bi 2 Sr 2 CaCu20 x (BSCCO), are composed of su- 
perconducting cuprate layers coupled by Josephson in- 
teraction. This system possesses the Josephson effects at 
the atomic scale ("Intrinsic Josephson Effect"). A rich 
spectrum of classical dc and ac Josephson phenomena 
have been observed in this system, see reviews 11121 

In a bulky superconductor the magnetic field applied 
along the layers generates a triangular lattice of Joseph- 
son vortices. The anisotropy factor 7 and the inter- 
layer periodicity s set the important field scale, B CI = 
<J>o/(2tt7s 2 ) (~ 0.5 tesla for BSCCO). When the mag- 
netic field exceeds B CT the Josephson vortices homoge- 
neously fill all layers. 3 Strong coupling between the vor- 
tex arrays in neighboring layers mediated by the in- 
plane supercurrents^ determines the static and dynamic 
properties of the lattice. Dynamic properties of the 
Josephson-vortex lattice in BSCCO have been exten- 
sively studied by several experimental groups (see, e.g., 
Refs. l5l6l7l8l) . 

When an external transport current flowing across the 
layers exceeds the critical current, the Josephson vor- 
tex lattice starts to move. In a homogeneous junction 
the critical current is determined by interaction with the 
boundaries. The simplest and most known case is a single 
small junction without inhomogeneities, where the field 
dependence of the critical current is given by the Fraun- 
hofer dependence, J c ($) = I c0 | sin(7T$/<l>o)|/(7r<i>/<i>o), 
with <I> being the magnetic flux through the junction. 
Observation of this dependence has been considered as 
an important confirmation of the dc Josephson effect. 9 
The same dependence is also expected for the junction 
stack with the lateral size smaller than the Josephson 
length. 10 In a single long junction the critical current has 
the rather complicated field dependence due to multiple 
coexistent states of the lattice.^ 

In the previous paper [12] we considered the behavior 
of the critical current for the dense Josephson-vortex lat- 



tice in a homogeneous wide stack for which the critical 
current is caused by interaction with the boundaries. We 
found that the boundary induces an alternating deforma- 
tion of the lattice. Averaging out the rapid phase oscilla- 
tions, we obtained that the lattice deformation obeys the 
sine-Gordon equation and decays inside superconductor 
at the typical length Lb/V&>, which is larger than the 
Josephson length, A,/ — 7s, and increases proportional 
to the magnetic field, Lb = XjB / B CI P^The stack is in 
the wide-stack regime if its lateral width L is larger than 
this typical length. In this situation the surface defor- 
mation and the total current flowing along the surface 
is uniquely determined by the lattice position far away 
from the boundaries. The surface current has oscillat- 
ing dependence on the lattice displacement and, due to 
the triangular-lattice ground state, the period of this de- 
pendence is half the lattice spacing. The total current 
flowing through the stack is given by the sum of two in- 
dependent surface currents flowing at the sample edges. 
The magnetic field determines the magnitude of the max- 
imum surface current (it is inversely proportional to the 
field) and sets the phase shift between the oscillating de- 
pendences of the two surface currents on the lattice posi- 
tion. One can trace that, due to the half-lattice-spacing 
periodicity of the surface current, a full period change of 
this phase shift corresponds to the change of the mag- 
netic flux through one junction, equal to the half flux 
quantum, <I>q/2. As a consequence, the maximum current 
through the stack has oscillating field dependence, which 
resembles the Fraunhofcr dependence: it has strong oscil- 
lations and overall 1/B dependence. However, the period 
of these oscillations is two times smaller: it corresponds 
to adding one flux quantum per two junctions and the 
critical current has local maxima at $ = fc<I>o/2. 

Oscillations of the flux-flow voltage in BSCCO mesas 
at slow lattice motion have been observed by Ooi et a?|5] 
The oscillations have the period of $o/2 per junctions 
and are caused the size-matching effect described in the 
previous paragraph. These oscillations have been re- 
produced by numerical simulations!^ More recently the 
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flux-flow oscillations in BSCCO mesas have been repro- 
duced and studied in more details by several experimen- 
tal groups.^- 5 * 16 * 17 * 18 ! Similar oscillations also have been 
observed in underdoped YBa2Ca306+a;P^ Size depen- 
dence of oscillations has been systematically studied in 
Refs. 115)171 It was found that at smaller lateral sizes 
and/or higher magnetic fields the crossover to the new 
oscillation regime takes place, in which the period be- 
comes <&o per junction, as in a single junction. 

Being motivated by recent experiments, in this paper 
we extend our consideration to the regime when the junc- 
tion size L is comparable with the length Lb and the 
system crosses over from the wide-stack to narrow-stack 
regime. As the length scale Lb increases with the mag- 
netic field, it also sets the field scale B L = B CZ L/Xj, at 
which this length becomes of the order of the junction 
length L. Therefore for a junction of a given size the 
crossover to the narrow-stack regime can be driven by the 
magnetic field, as it was observed experimentally.^ i>o/2- 
periodicity of the critical-current oscillations holds un- 
til interactions between surface deformations can be ne- 
glected. This interaction becomes progressively stronger 
with decreasing the ratio L/Lb- Surface deformations 
at L > Lb can be described as partial sine-Gordon 
solitonsPS The relative sign of two solitons at the op- 
posite edges is determined by the magnetic field and 
the lattice positions. At the integer-flux-quanta points 
$ = k$ the surface solitons always have the same 
sign and repel each other. As a consequence, the am- 
plitude of surface deformations drops and the critical 
current decreases. At the half-integer-flux-quanta points 
$ = (fc + 1 /2)<& situation is the opposite: the surface 
solitons always have opposite signs and attract each other 
leading to enhancement of the surface deformations and 
increase of the critical current. Therefore the interaction 
between the surface solitons leads to the crossover be- 
tween the <&o /2-periodic oscillations of the critical current 
and <i>o-periodic oscillations. This crossover occurs via 
suppression of the current peaks at the points $ = fc$o 
and enhancement of the current peaks at the points 
<3> = (fc + l/2)$o- Such behavior is consistent with recent 
studies of the oscillatio ns of the flux-flow voltage in mesas 
with small lateral sizes IMliSUZI The crossover in t he vol t- 
age oscillations also has been studied numerically ! 16 * 2 1 I 

In the region L ~ Lb the lattice structure is deter- 
mined by competition between two energies: the interac- 
tion with boundaries and the bulk shearing interaction 
between the Josephson- vortex planar arrays in neighbor- 
ing layers. The interaction with the boundaries favors the 
aligned rectangular arrangement of the Josephson vor- 
tices while the local shearing interaction favors the tri- 
angular lattice. The boundary interactions decay slower 
with increasing field then the shearing interaction and 
become dominating at large fields. On the other hand, 
the boundary interaction energy has oscillating field de- 
pendence and vanishes at the points $ = k$Q. At these 
points the shearing interaction is relevant at any mag- 
netic field. In addition, the interaction with the bound- 




FIG. 1: Left plot: Phase diagram of the Josephson-junction 
stack in the coordinates [reduced junction size Lj A j]- [reduced 
field h — B / B cr \. The solid lines mark fields corresponding to 
the integer values of the magnetic flux per junctions, $ = fc$o- 
Dotted lines show boundaries of the transitions into the rect- 
angular lattice for ground state. Right plot: The same dia- 
gram in the coordinates [number of flux quanta per junction 
$/$o]- [ratio L/L B = L/(hXj)]. 



aries is suppressed by the external current. 

In the region L -C Lb the rectangular arrangement of 
vortices is realized in most part of the phase space and 
the field dependence of the critical current approaches 
the classical Fraunhofcr dependence. Two important de- 
viations persist at all fields and sizes: (i) Near the points 
$ = fc^o the phase transitions to the triangular lattice 
always take place. Critical current at these points never 
drops to zero and actually always has a small local max- 
imum; (ii) Away from the points $ = fc<I>o the critical 
current is reached at the instability point of the rectan- 
gular vortex lattice and it is always somewhat smaller 
than the "Fraunhofer" value, 7 c0 | sin(7r4>/$ )|/K<iV$o|- 
Therefore, we somewhat revise the statement of Ref. [TD] 
that behavior of a small stack is identical to a single small 
junction. The described results has been summarized in 
a short paper published in the conference proceedings.^ 
The static lattice structures at different fields and sizes 
have been also studied numerically in Refs. [23] and [5T] 
and the results are in a qualitative agreement with the 
described picture. 

To illustrate a general picture, we show in Fig. [l] (left) 
the ground-state phase diagram of the junction stack in 
the field-size plane, where h = B / B cr is the reduced mag- 
netic field. Solid lines mark the magnetic fields corre- 
sponding to integer values of the magnetic flux per junc- 
tion and dotted lines show boundaries of the rectangular- 
lattice regions for zero current through the stack. As one 
can see, these regions appear above the line h — 1.48L/Xj 
and they are always bounded by the integer-flux quanta 
lines. In the right plot the same diagram is replotted 
in different coordinates, number of flux quanta per junc- 
tion Q/$o vs ratio L/Lb = L/(hXj), which controls the 
wide-stack/narrow-stack crossover. It is interesting to 
note that in these coordinates the diagram is periodic 
with respect to the magnetic flux through the junction. 

The paper is organized as follows. In Section [TT] we 
outline derivation of the phase distribution in the case 
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of alternating from layer to layer solution. Averaging 
with respect to the rapidly changing phase oscillations, 
we obtain equation and boundary conditions for the slow 
lattice deformation. We also express the lattice energy 
and current flowing through the stack via this deforma- 
tion. In Section [llT| we obtain and analyze solution for the 
smooth phase in terms of elliptic integrals. We found that 
the problem can be reduced to solution of three nonlin- 
ear coupled equations for three unknowns, the boundary 
phases and elliptic parameter. In Section |IV| we derive 
a criterion for the transition into the rectangular-lattice 
state. In Appendix [B] we consider weak finite-size ef- 
fects in the wide-stack regime and analytically compute 
exponentially small finite-size corrections to the critical 
current, which break the <j> /2-periodicity of oscillations. 
In Section [V] we present results of numerical analysis of 
the crossover between the wide-stack and narrow-stack 
regimes with the increasing magnetic field. We obtain the 
oscillation patterns of the critical current for stacks with 
different lateral sizes and find location of the rectangular- 
lattice regions in the current-field plane. In Section VI 
we reanalyze in detail the narrow-stack regime using in- 
dependent analytical approach. In Section |VII| we con- 
sider the voltage oscillations in the case of slowly mov- 
ing lattice and relate these oscillations with the critical- 
current oscillations. We elaborate the recipe to extract 
the anisotropy factor from the voltage oscillations. 



II. PHASE DISTRIBUTION AND ENERGY OF 
FINITE STACK ASSUMING ALTERNATING 
SOLUTION 

We consider a Josephson-junction stack consisting of N 
layers, N 3> 1, with lateral size equal to L in a magnetic 
field B > B cr applied along the layers. At high magnetic 
fields one can neglect screening effects. In this case the 
stack is described by the energy functional (per layer and 
per unit length in the field direction) of the layer phases 
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where J is the in-plane phase stiffness and Ej is the 
Josephson energy per unit area. To simplify deriva- 
tions, we introduce reduced coordinate, u = x/Xj, with 
Xj = \J J/Ej, reduced magnetic field, h = 2ttsXjB/$ , 
and reduced energy £ = E/(EjXj). We also represent 
the phase variable in the form, which naturally describes 
the dense triangular lattice in the bulk in the limits h 1 
and L/Xj ^> h, (p n {u) = 4> n (u) + an + irn(n—l)/2, where 
the phases 4> n {u) are assumed to be small and rapidly 
oscillating and the parameter a will describe lattice dis- 
placement. The the reduced energy per one junction and 



per unit length in the field direction can now be repre- 
sented as 



£[</>n] =JjY1 / dU 



1 



2 V du 

— cos {4> n +i — 4> n — hu + a + 7m)] . (2) 

The oscillating behavior is determined by the reduced pa- 
rameter hL which is directly related to the total magnetic 
flux through one junction <f> = BLs, 

hL = 27T$/$ . 

We consider the stack containing a very large number 
of junctions, N ^> 1. This will allow us to focus on 
bulk behavior and neglect c-axis boundary effects com- 
ing from the top and bottom junctions, which give 1/N 
corrections to the bulk results. We will also not consider 
potentially interesting "parity effects" , small differences 
between stacks containing odd and even number of junc- 
tions, which have the same order. Our results are also not 
influenced by possible perturbations of the current distri- 
bution near the boundaries. Due to the large anisotropy 
of the material, the "bulk" current distribution usually 
is realized already in the second junction in the stack. 

Following Ref . [T2 we will assume the alternating phase 
distribution in the form <j) n = (— l) n <fi. This distribu- 
tion describes both the deformed triangular lattice in the 
wide-stack regime and the transition to the rectangular 
lattice in narrow-stack regime. Substituting this presen- 
tation into Eq. we represent the energy functional 
as 



£(a; 



)= ( du 
Jo 



du 



— sin (20) sin (— hu+a) 



(3) 



The lattice energy as function of displacement, £(a), is 
determined by the minimum of the functional £ (a; 4>{u)) 
with respect to </>(«), £{oi) — min^ [£{a; (f>)]. As the en- 
ergy functional has a symmetry property a — > a + n, 
<p — > —0, the energy £{a) is 7r-periodic with respect to 
shift of a, £ (a + it) — £(a). 

The ground-state phase distribution <fi obeys the fol- 
lowing equation 



du 2 



+ 2 cos (24>) sin (-hu + a) = 0, 



(4) 



which has to be solved with the boundary conditions 
d(t> 



du 



0, for u = 0, L. 



(5) 



In the limit /i » 1 further significant simplification is 
possible: we can average out rapid phase oscillations and 
derive a simplified equation for the smooth phase per- 
turbation. We split the total phase into the smooth and 
rapidly-oscillating components 



(f)(u) = v(u) + (j)(u), 



(G) 



4 



where we assume \<f>\, \dv/du\ <C 1. The maximum value 
of v(u), Umax = corresponds to rectangular lat- 

tice. The rapidly-oscillating phase by definition obeys 
the equation 



du 2 



+ 2 cos (2v) sin (— hu + a) = 0, 



which has the following approximate solution 



„ cos (2v) sin (— /m + a) . 



(7) 



(8) 



To the first order with respect to <j>, equation for v(u) is 
given by 



du 2 



4(f) sin (2v) sin (—hu + a) = 0. 



Substituting into this equation the oscillating phase ([8]) 
and averaging it with respect to rapid oscillations, we 
finally obtain the sine-Gordon equation for the smooth 
phase 



cPv 
du 2 



h 2 



sin (Av) = 0. 



(9) 



Computing the derivatives of the rapid phase ([8]) at the 
edges, we also derive the boundary conditions for the 
smooth phase, 

~(0) = % cos (2v a ) cos (a) , (10a) 
du h 

^-(L) = y cos (2vl) cos (—hL + a) (10b) 
du h 

with wo = v(0) and vl = v(L). The local current den- 
sity, j(u) = j j sm[6 n n+ i(u)] with jj being the maximum 
Josephson current density, is determined by the gauge- 
invariant phase difference, 6'„ i „+i(u) = 4> n +\ — 4>n — hu + 
a + 7m, which is related to v(u) as 



9 ra ,n+i ~ —hu + a + 7m — (—l) n 2v 



h 2 



cos (2v) sin (—hu + a + 7m) . 



(11) 



Substituting the phase presentation ^ and ^ into 
the energy ^ and averaging with respect to the rapid 
oscillations, we derive the energy functional in terms of 
the smooth phase v, 



E(pt; v) 



[sin (2vq) cos (a) —sin (2vl) cos (—hL + a)] 



du 



1 + cos 4w 
2h 2 



(12) 



To shorten notations, we omitted the arguments h and 
L in £(a,h, L;v). Minimization of this energy func- 
tional with respect to v(u) gives the energy as a func- 
tion of the lattice shift a, £(ot). Minimum of this energy 
with respect to a gives the ground state for given h and 



L. Higher-energy states at other values of a typically 
carry a finite current. The total Josephson current flow- 
ing through the stack is proportional to d£/da. Taking 
derivative of the functional (12 1 with respect to a, as- 



suming that at every a it is minimized with respect to 
v(u), we obtain 

J(a) = — [— sin(2i'o) sina+sin(2vi) sin(— hL+a)] . (13) 

The unit of current in this equation is jjXjw where jj is 
the Josephson-current density and w is the junction size 
in the field direction. An important consequence of this 
equation is that nonzero current exists only if the surface 
deformations vq and vl are finite. Further analysis is 
based on Eqs. @, ([lOj) , @, and 



III. SOLUTIONS FOR SMOOTH PHASE 

A general solution of the sine-Gordon equation ^ can 
be found in terms of elliptic integrals. From the first 
integral of Eq. ^ we obtain 



dv /-JXjm — cos 2 (2v) 

— = d d v2 

du h 



(14) 



with 5d = sign[efo/du] = ±1 and m is the elliptic parame- 
ter which has to be found from the boundary conditions. 
From this equation we obtain implicit equation for defor- 
mation v(u), 



dv 



v a \J\jm — cos 2 (2v) 



= 6 d V2u/h. (15) 



To rewrite this equation using standard elliptic integrals, 
we introduce a new variable if, 



tp = 7r/2 + 2v. 



(16) 



This variable has its own physical meaning: it describes 
the alternating deformation of the interlayer phase dif- 
ference with respect to the rectangular-lattice state. In 
particular, ip = corresponds to the rectangular lattice. 
Using these variables, we can rewrite Eq. (fl5| as 



where 



m[F(<p,m)-F(<p ,m)] = 5 d V8L/h (17) 



F(tp,m) 



di 



o \J 1 — m sin 2 x 



is the incomplete elliptic integral of the first kind.^ 

In the limit of very large L the deformation v has to 
vanish far away from edges meaning that m — > 1. For 
finite-size junctions, depending on hL and a, the solu- 
tion v(u) can be either monotonic or nonmonotonic. For 
the nonmonotonic solution the derivative dv/du and the 
parameter 5 d change sign inside. The monotonic solu- 
tion can either change sign (vqVl < 0, m < 1) or not 
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(voVl > 0). For large L the monotonic solution corre- 
sponds to the two surface partial solitons of the same 
sign and the nonmonotonic case corresponds to the sur- 
face solitons of opposite signs. A mathematical structure 
of solutions for these two cases is different and we will 
consider them separately. 

First, we find some general relations between the 
boundary phases and the parameter m. From the bound- 



ary conditions ( 10a I, ( 10b I , and Eq. ( 14 1 we obtain equa- 
tions 



5 \/l/m — cos 2 (2v ) = \/2cos (2v ) cos (a) , (18a) 
S L yjl /m-cos 2 (2v ) = V2 cos(2v L )cos(hL-a) (18b) 

with 6q = 6d(0) and Sl — $d(L) (for monotonic solution 
5q = Sl and for nonmonotonic solution So = —Sl)- As 
\ v 0,l\ < 7r /4, the inequality cos(2u 0j l) > always holds 
meaning that <5 = sign[coso;] and = sign [cos (hL — a)]. 
From the conditions (18 1 we obtain the following results 



for the boundary deformations 



cos (2«o) = 



cos (2vl) 



1 1/m 
1 + 2 cos 2 (a) 



1/m 



1 + 2 cos 2 (hL - a) 



(19a) 
(19b) 



or, in terms of variable tp (16) 



sin tpo = 



an (fi L = 



1/m 



1 + 2 cos 2 (a)' 
1/m 

l + 2cos 2 (hL-a)' 



(20a) 
(20b) 



From these considerations one can conclude that the 
monotonic solution realizes for all cc's for hL = 2kn (mag- 
netic flux per junction equals to integer number of flux 
quanta, $ = fc^o) and the nonmonotonic solution real- 
izes for hL = (2k + l)ir (magnetic flux equals to half- 
integer number of flux quanta, $ = (k + 1/2) $ ). If the 
magnetic flux through one junction is not close to these 
special values then the solution changes from monotonic 
to nonmonotonic depending on the lattice phase shift a. 
Location of different types of solution depending on hL 
and a is illustrated in Fig. [2] 

One can distinguish two important special cases cor- 
responding to symmetric locations of the lattice with re- 
spect to the boundaries. The first case, a = hL/2 + irk. 
vl = —Vo, corresponds to the monotonic solution and 
the second case, a = hL/2 + n /2 + irk, vl — vq, corre- 
sponds to the nonmonotonic solution. The energy al- 
ways has extremums at these values of a. Moreover, a 
detailed study shows that the ground state is always re- 
alized at one of these values of a. At large L the system 
switches between these locations in the vicinity of the 
points hL = (2k + l/2)n, as it is illustrated in Fig. [2] 




22kn hL=2:iO/cD n (2k+1)7r 



u 



FIG. 2: Regions of different types of solution for the smooth 
alternating deformation v(u) in the plane hL-a for L/h ~ 
1. Typical solutions are illustrated for five marked points. 
Dark grey color marks regions of monotonic solution changing 
sign (3), light grey color marks regions of monotonic solution 
not changing sign (2 and 4), and white regions correspond to 
nonmonotonic solution (1 and 5). Light-grey regions shrink 
with increasing L. The black line illustrates the location of 
the ground state. 



In the vicinity of these switching points the energy has 
minima at the both values of a. 



In general, the conditions ( 19a I and ( 19b I are not suffi 



cient to determine signs of the edge deformations and 
vl- In the limit of large L the deformation v(u) has to de- 
cay from the edges leading to relations sign[uo] = — So — 
— sign[cosa] and sign[t>i] = Sl = sign[cos (hL — a)]. In 
this case we also obtain conditions 



tan(2v ) = -So^/m (2 cos 2 (a) + 1) - 1, 
tan(2w L ) = S L \Jm (2cos 2 (hL-a) + 1) - 1 



which fix signs of vq and Vl- For large values of L/h 
monotonic solution typically changes sign inside. How- 
ever for finite L there are intermediate regions exist lo- 
cated near lines a — n/2 + nk and a = hL + n /2 + irk 
where solution is still monotonic but does not change 
sign. We now proceed with analyzing separately the 
monotonic and nonmonotonic solutions. 
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A. Monotonic solution 

The monotonic solutions realize for ranges of a where 
cos a cos (hL — a) > (grey regions in Fig. [2|. For such 
solution we obtain from Eq. ( 15 I relation connecting the 



parameter to with the boundary deformations 



C 

J Vn 



dv 



a/1/to — cos 2 (2v) 



sign[cosa]V2L//i. (21) 



Using previously introduced variable p (16 1, we can 
rewrite Eq. pTj) via elliptic integrals as 



to [F(ip L , to) — F(ip , to)] = sign[cosa]\/8i//i (22) 



This equation together with boundary conditions (20a I, 



and (20b) have to be solved to find the three unknown 



constants (po, <pl, and to, which completely determine the 
solution. The boundary deformations vq = (po~ 7r /2)/2 
and Vl = ( l Pl~' k /2)/2 may have either the same sign or 
opposite signs. In Appendix[X]we find the boundary sep- 
arating these two types of solution (boundaries between 
dark-grey regions and light-grey regions in Fig. 21). 

The energy ( 12 1 and Josephson current (113 ) can be 
represented as 



£ = — — [cos (ip ) cos (a) — cos (ipl) cos (hL — a)] 



1 



I2mh 



\E(ip L ,m) - E(cp ,m)\ 



mh 2 



(23) 



J(a) — \ [cos (ipo) sin a — cos (ip^) sin (—hL + a)] , (24) 
h 



where E(tp,m) = f? yl — msiv^x dx is the incomplete 
elliptic integral of the second kindP^l 



Adding these two equations, we obtain equation connect- 
ing v m with the boundary deformations vq and vl 



dv 

a/cos 2 (2v m ) — cos 2 (2v) 
' m dv 

yjcos 2 (2v m ) — cos 2 (2?;) 



S d V2L/h. (25) 



This equation together with boundary conditions ( |19a l 
and (19b I, represents the full system for determination 
of three unknown constants vq, vl, and to= 1/ cos 2 (2w m ). 
To rewrite these equations in terms of elliptic functions, 



we again transfer to variable (16 1. Then equation (25 1 
can be rewritten as 

V&L 

\fm [2F((p m , m)-F(tp Q , m)-F((p L ,m)] = 5 d —^—. (26) 

with ip m = tt/2 + 2v m . The boundary conditions in 
terms of ipo and (p^ are again given by Eqs. (20a I 



and ( 20b I . The elliptic parameter is related to (p m 
as to = 1/ sin 2 tp m leading to the following relation, 
F(tp m ,m) — K(l/m) j \/rri. The energy for nonmono- 
tonic solution can be represented as 

£ = — — [cos ipo cos (a) — cos tp L cos (—hL + a)] 



1 



'2mh 



\2E(cp m ,m)-E((p ,m)-E(ip L ,m) | - 



ih 2 



(27) 



and the Josephson current is again given by Eq. (24 1 



One can easily check that the nonmonotonic solution 
matches the monotonic solution at the boundaries. For 
example, for cos a = the extremum is located at the 
boundary u = 0. In this case we have 4>o = <p m = 



arcsin(l/y / ro) (or 7r — arcsin(l/ v /ro)) and Eq. (26 1 co- 
incides with Eq. ([22 1. 



B. Nonmonotonic solution 

The solution is nonmonotonic in the regions given by 
cos a cos (hL — a) < (white regions in Fig. [2]). In this 
case the function v(u) has extremum at some point u — 
u m so that we can rewrite Eq. ( 14 1 as 



dv / - y /\/m-cos 2 (2v) 

— = S d V2^ , for u < u m , 

du h 

dv r \l\jm — cos 2 (2v) 

— = -5 d V2 , for u m <u < L 

du h 

with S d = sign [cos a]. As dv/du = at u = u m , we have 
to= 1/ cos 2 (2i> m ) > 1 with v m = v(u m ). Integrating the 
first equation from to u m and the second equation from 
L to u m , we obtain 



dv 



a/1 /to — cos 2 (2v) 
dv 

\J\jm — cos 2 (2v) 



= S d V2u m /h, 

= S d V2 (L - u m ) jh. 



C. Alternative presentation of equations via 
elliptic integrals 

Using known relations for the elliptic integrals 
F(tp, m) = —^=F(ijj, 1/to) with sin ip = a/to sin ip, (28a) 

V TO 



E(ip, m) = — = [(1 — m)F (tp, 1 /m)+mE (ip, 1/to)] 

a/To 



(28b) 

valid for tp < n/2 and sin <p> < 1/a/to, one can rewrite 
equations (p2| and ( 20 1 for the case of the same-sign 



monotonic solution in the following alternative form 

F (ip L ,fh) — F (ipo,m) — sign[cos a]V$>L/h, (29a) 
1 



sin ipo = 
sinip L = 



a/1 + 2 cos 2 (a)' 
1 

a/1 + 2 cos 2 (hi- a/ 



(29b) 
(29c) 
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with m = 1/m. Correspondingly, the energy (23 1 in this 
representation is given by 



£-- 



1 — to sin i/>ocosa— y 1 — to sin ipLCOs(hL — a) 



' \E(ip ,fh)-E(ip L ,m)\- (2 ™ )L . (30) 



V2h 



In the case of nonmonotonic solution, using relations 
(28a I and (28b I for the elliptic integrals, one can rewrite 
equation ( |26[ ) in the following equivalent form 

2K (to) - F {ip , to) - F (il> L ,m) = S d V&L/h. (31) 



where to = 1/m, ip and tpL are given by Eqs. ( |29b[ ) and 
( 29c ) , and represent the energy as 



\J 1 — to sin 2 V>o |cos a| 
+ \/ 1 — to sin 2 |cos (/iL — a) | 



// 2 



(32) 



This representation is especially useful in the case of 
large to (small in). In particular, it will allow us to study 
transition to the rectangular-lattice state corresponding 
to the limit to — > oo, which we will consider in the next 
section. 



IV. TRANSITION TO THE RECTANGULAR 
LATTICE 

An important particular case is the solution of Eq. |9| 
corresponding to rectangular lattice, v — ±7r/4 or Lp = 
0, 7r. This case corresponds to the limit to — > oo (to — > 0). 
The energy of the rectangular lattice coincides with the 
well-known result for a single junction 



£rect (a) 



— sin | a 
2 J 



hL 
~2 



(33) 



and has minimum £ rcc t — — 2 |sin(/ii/2)| /h at a = 
hL/2 + 5n/2 with S = sign [sin (hL/2)}. 

To find condition for the transition to the rectangular 
lattice, we take the limit to — > in Eq. (31 1 for non- 
monotonic solution. Using relations K(0) = tt/2 and 
F(ip, 0) = -0, we obtain 



7T - 1p - 1p L = V&L/h, 



(34) 



where ipo an d ipL are given by Eqs. (29b I and (29c I. Us 



ing these definitions, the condition for the rectangular 
lattice can be rewritten in an explicit form as 



\/2 (j cos (hL — a) \ + \ cosa|) 
Vd^cos 2 a) (l + 2cos 2 (hL-a)) 



<sm 



V8L N 



(35) 



o=ko r 
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FIG. 3: Typical lattice deformations at large L for $ = fc$o 
(upper plot) and $ = (k + l/2)"I>o (lower plot). In the former 
case the surface partial solitons have the same sign and repel 
each other while in the latter case they have the opposite signs 
and attract each other. The insets in both plots illustrate cor- 
responding displacement fields of the Josephson-vortex lattice 
in the two neighboring layers. 



This equation gives the transition criterion in general 
case, including the current-carrying states. In particular, 
the rectangular lattice gives a local energy minimum at 
a = hL/2 + ir/2 in the regions \hL/2ir- (k + 1/2) | < 1/4 
if the inequality 



\sin(hL/2)\ < tan (V2L/h 



(36) 



is satisfied. The rectangular lattice first appears in the 
ground state at points hL = (k + l/2)27r for L/h <h = 
arctan (V2) /V2 ~ 0.675. This value is marked in the 
right plot of Fig. [T] 



WIDE-STACK/NARROW-STACK 
CROSSOVER 



In this section we investigate in detail the crossover 
between the wide-stack and narrow-stack regimes. As 
this crossover is driven by the reduced parameter h/L, 
for a junction with size L the crossover takes place with 
increasing magnetic field at size-dependent field Bl = 
L$ /(27r 7 V). 

At large L, L h or B <C -Bl, the smooth alternat- 
ing deformation has solutions in the form of two isolated 
surface solitons l 2 ^ The monotonic solution corresponds 
to the solitons of the same sign and the nonmonotonic 
solution corresponds to the solitons of opposite signs, as 
it is illustrated in Fig. [3] If one neglects the interaction 
between the solitons then the relative sign of surface soli- 
ton has no importance and the total Josephson current 
is given by the sum of two independent surface currents, 
which do not depend on the soliton signs As a con- 
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FIG. 4: The representative dependences of the energy (upper 
left plot) and Josephson current (lower left plot) on parameter 
a for L — 4 and three values of flux per junction = 3.1, 

3.25, and 3.4 within one oscillation period (the curves are 
marked by these values) corresponding to the magnetic fields 
h = 4.87, 5.105, and 5.34. The parameters correspond to the 
crossover region L/ft~ 1. Dark grey marks region of mono- 
tonic solution and light grey marks region of nonmonotonic 
solution. Arrows in the energy plot mark values of a corre- 
sponding to ground state. Pictures in the right column illus- 
trate structure of ground states at <&/<&o = 3.1 and 3.4. The 
lines show oscillating z-axis currents in neighboring layers and 
ellipses mark centers of the Josephson vortices. 



sequence, the product hJ c has periodicity of half flux- 
quantum per junction and reaches maxima at hL = nj 
($ = j$ /2) with hJ c ps 1.035. At finite L the interac- 
tion between the surface solitons disturbs such period- 
icity. At large L one can derive analytically corrections 
to the infinite-I/ results, see Appendix [B] for details. In 
particular, near the maxima hL = Ttj ($ = j<&o/2), we 
find the finite-size correction, 



5J c (h, nj) 



4.544 



(-l)j 



h 



exp 



V8L 



(37) 



As we can see, the finite-size effects increase the critical 
current maxima at 4>= (fc+l/2)<I>o (j=2fc+l) and reduce 
the critical current maxima at <I> = fc<I>o (j = 2fc). In the 
wide-stack regime, however, these corrections are expo- 
nentially small, which explains nice $o/2-periodic oscil- 
lation of the flux-flow voltage observed in this regimeP^ 
In the whole range of fields and sizes we explore the 
phase diagram numerically. To find the ground state and 
the critical current at given h and L, we study depen- 
dences of the lattice structure, energy, and Josephson 
current on the lattice phase shift a. First, we have to 
find the boundary deformations <po, fL and the ellip- 
tic parameter m using the boundary conditions ( 20 1 to- 



gether with either Eq. (22 1 for cos a cos ( h L — a) > or 
Eq. (26 1 for cos acos(hL— a) <0. Using obtained values, 
we compute the energy from Eq. ( 23 1 or ( 27 1 and the 
current from Eq. ( 24 1 . This procedure has been imple- 



mented in Mathematica. Figure [3] shows representative 
a-dependences of the energy and current for L = 4 and 
three values of flux per junction 4>/<&o = 3.1, 3.25, and 




0123456 2 4,6 

° /<J5 ° h/L 

2 0.5 1 1.5 2 




FIG. 5: The field dependences of the critical current for three 
different sizes L = 3, 4, and 6. To emphasize the periodic 
nature of dependences, we plot the product $ J c (in units of 
Jj$o/27T where Jj is the total maximum Josephson current) 
vs &/&o- One can observe the crossover from <E>o/2-periodic 
oscillations to "l>o-periodic oscillations with increasing $. At 
larger L the crossover takes place at larger $/$o- The dashed 
line in the L = 3Aj plot shows function 2| sin(7r$/$o)| corre- 
sponding to usual Fraunhofer dependence in small Josephson 
junctions. Top axes show the parameter h/L (in real units 
h/L = B/Bl with B L = L$ /(27r7 2 s 3 )). The crossover al- 
ways takes place at the same value of ratio h/L ~ 1.5. In 
the plots for L = 4Aj and 6Aj dotted and dashed lines show 
universal dependences of product $ times current maxima at 
$ = (k + l/2)$ and $ = fc<I>o ($>J max ,2) on the 

parameter h/L. 



3.4 within one oscillation period. The minimum of the 
energy with respect to a determines the ground state 
and the maximum of the current determines the critical 
current. 

Figure [5] shows the field dependences of the critical 
current for three different values of the junction size L, 3, 
4, and 6. One can observe that with increasing field $0/2- 
periodic oscillations smoothly transform into <i>o-periodic 
oscillations. This occurs via suppression of the peaks at 
$ = £;<i>o and enhancement of the peaks at $ = (fc + 
l/2)$o- Such behavior of the critical current has been 
recently reproduced by numerical simulations by Iric and 
OyaP3 Experimentally, the crossover between Qq/2- and 
<I>crPeriodic oscillations has been observed in the flux- 
flow voltage by Kakeya et ali^ The crossover field can 
be arbitrarily defined as a field at which a fc<i>o-peak in 
the product <!>J C drops below the half of a (fc + 1/2)<J> - 
peak. At larger L the crossover takes place at larger 
field and larger ^/^o but it always occurs at the same 
ratio h/L, h/L ~ 1.6. Important property of the system, 
discussed in Sec. |IV| is the transformation of the lattice 
into the rectangular state at sufficiently large h/L. This 
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FIG. 6: The field dependences of the critical current (solid 
lines) for sizes L = 2.5, and 4 for the same range of the ratio 
h/L, h/L < 3.5 shown on the top axes. Shaded areas show the 
regions of stable rectangular lattice. The rectangular-lattice 
regions first appear in the vicinity of points $ = (k + 1/2) $o 
when h/L exceeds l/h w 1.48. When h/L exceeds l/h = 
2.07 the rectangular lattice remains stable at these points up 
to the critical current. 



property is illustrated in Fig. [6] where behavior of the 
critical current is shown together with regions of stable 
rectangular lattice in the field-current diagram for two 
sizes, L = 2.5 and 4. 

Lets consider in more detail behavior of the critical- 
current maxima at $ = (j/2)<J>o- We start with the case 
of half-integer flux quanta per junction, <f> = (k + l/2)$o 
(hL = (2k + 1) 7t). In this case a nonmonotonic solution 
is realized for all a. As for sufficiently high h/L the lattice 
transforms into rectangular state, it will be convenient 
to use presentation given by Eqs. 



(29b) 



and (31) 



which naturally describes this transition. 
Eq. (34), in the range 



cos a > tan (V2L/h) /V2 



is follows from 



(38) 



the rectangular lattice is realized for which fh — and 
the energy and current are given by 

2 

£i(a) = — — |cos a I 
2 

Ji(a) = — sin a sign[cosa] 

In the opposite range, cosa < tan. (y2L/ti\ /y/2, the 
solution has the form of deformed lattice. In this 
case we have Vo = i>L = arcsin(l/-\/l + 2cos 2 a) (or 



7T — arcsin(l/Vl + 2cos 2 a) ) and Eqs. (31 1 becomes 

1 



K (fh) — F arcsin 



\/l + 2 



cos^ a 



, to 



'2L 



(39) 



Solving this equation with respect to the elliptic param- 
eter fh(a), we can obtain the energy and current 



£i(a) 



V2 



+ ^-[E(m)-E(^ 0l m)] 
h 



in 

^--^ cos (a) 

1 + 2 cos 2 (a) v ; 

(2-m)L 



h 2 



(40) 



2 / fh 

Ji(a) = ^v 1 "^^^) sin(a)sign[cos(a)] (41) 



The critical current at $ = (k +l/2)<f>o is given by 
>^max,i — max Q [J 1 (a)]. At small L, L/h <C 1, the max- 
imum critical current is realized at the instability point 
of the rectangular lattice cos (a) ~ L/h giving 



J 11 



h 



2h 2 



(42) 



It is always somewhat smaller than the "Fraunhofer" 
value 2/h. 

It was obtained in Section |IV| that the rectangular 
lattice is realized in ground state (a — 0) at points 
$ = (k + l/2)$ for L/h < h = 0.675. If, however, 
L/h is only slightly smaller than this value, the rectan- 
gular lattice becomes unstable with increasing current 
and the configuration at the critical current still corre- 
sponds to the deformed lattice. We found that there is 
another typical value of the ratio L/h, L/h = I2 ~ 0.484, 
below which the rectangular lattice remains stable up to 
the critical current. Both typical values of h/L, l/l\ and 
I//2, are marked in Fig. [6] One can see in that for both 
shown stack sizes, L = 2.5 and 4, the rectangular lat- 
tice first appears around points $ = (k + 1 /2)<&o when 
h/L exceeds l/l\ and its stability range extends up to 
the critical current when h/L exceeds l/h- 

For integer flux quanta, <f> = k$ (hL — 2irk), a 
changing-sign monotonic solution always realizes, vl — 
—vq. In the case cos a > this corresponds to cpn = 
7r — tpL — arcsin[m(l + 2cos 2 (a))] -1 / 2 and Eq. ( 22 1 can 
be reduced to the form 



m (K (to) - F ((j> , to)) = V2L/h. 



(43) 



Solving this equation with respect to to, we can obtain 
the energy from Eq. (|23[) and current from Eq. ( 24 1 



2 V^2 L 

£2 = -- cos (j) cos a+ [E(m) -E((f> ,m)]- 



J2 (a) = — cos 4>o sin a 



(44) 
(45) 
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The critical current at <f> = (fc + l/2)<I>o is given by 
>/max, 2 = max Q [J2(a)]. At small L inequality cos^o <C 1 
holds. In this limit, using relation F(<pQ,m) ~ K(m) — 
(n/2 — (f>o) / \J\ — m, we can approximately rewrite Eq. 
(431 as 7r/2-0o = ( \f2L/h) x /\~\/m. As sin^o is 



close to one, Eq. (20a I gives 1/m ~ 1 + 2 cos 2 (a) and 
o = 7r/2 — (2L/h) cos a. Therefore we obtain for the a- 



dependent current ( 45 1 



, . 2L . n 
Jzipt) ~ 717 sin 2a. 



(46) 



The maximum is realized at a — n/A giving the following 
result for the critical current 



2L/h 2 



(47) 



i.e., it decays at large h as l/h 2 but never drops to zero 
as for usual Fraunhofer dependence. The behavior in the 
narrow-stack regime will be considered in more details in 
the next section. One can see that the critical currents 
at both maxima J max ,a (a = 1,2) have the same scaling 
property: the product hJ maxa depends only on the ratio 
L/h. These scaling dependences are plotted in Fig. [5] in 
the plots for L — 4 and 6. 



From this energy we obtain equation for f(u), 

^ + -sm(2^) = 0, 

and the boundary conditions 

df . . 4 , / „ hL 



(49) 

(50a) 
(50b) 



For small L Eq. ( 49 1 can be solved as expansion with 



respect to powers of u — L/2, 



f = <Pa + alu- 



2(u-L/2Y 
h 2 



sin(2^ a ) (51) 



with ip a = f(L/2). Boundary conditions (50 1 give two 
equations for two unknown variables, the midpoint phase 
ip a and the linear slope a. We obtain two types of 
solutions: (i) the rectangular-lattice solution a = 0, 
sin ip a = and (ii) the deformed-lattice solution. In the 
leading order with respect to the small parameter L/h, 
the latter solution can be represented as 



VI. NARROW-STACK REGIME 

Lets consider in more detail the narrow-stack regime at 
L/h <C 1 . In this regime interaction with the boundaries 
is typically stronger than the bulk shearing interaction. 
As a consequence, the boundaries stabilize the rectangu- 
lar lattice configuration in most part of the phase dia- 
gram. The exception is the narrow regions in the vicin- 
ity of the integer-flux-quanta points <& = fc<&o where the 
interaction with the boundaries vanishes and the rectan- 
gular lattice looses its stability. The rectangular lattice 
also becomes unstable near the critical current. In this 
section we will study in details this behavior. Instead of 
using asymptotic behavior of elliptic integrals, it is more 
transparent to use as a starting point the equation for 
smooth alternating deformation the boundary con- 
ditions (10) and the energy (12 1. It will be more con- 



venient to use variable if given by Eq. ( 16 ) (instead of 
v) from the very beginning, because it vanishes in the 
rectangular-lattice state. We also introduce a new vari- 
able for the lattice phase shift, 

[3 = a + Ti/2-hL/2, 

which will facilitate a more compact presentation of re- 
sults. In terms of the variables ip(u) and p the energy 
( 12 1 can be rewritten as 



. , , hL\ / hL\ 

cos ip sin [ p+ — J -cos tp L sm I p- — 1 



du 



dip 
du 



1 — cos(2tp) 
2h 2 



(48) 



a » — sin (if a ) sin /3 cos (hL/2) , (52) 



COS (fa) 



h sin(hL/2) cos j3 
Ll + 2sin 2 Pcos 2 (hL/2)' 



(53) 



As follows from the last equation, the deformed-lattice 
solution does not exist if 



h | sin (hL/2) cos p\ 
Ll + 2sin 2 P cos 2 (hL/2) 



> 1. 



(54) 



In this case the configuration must be the rectangular 
lattice. The solution (53 1 also includes the case of the 
ideal triangular lattice f a = Tr/2 which is always realized 
if either sin (hL/2) = or cos/3 = 0. As we consider the 
region L/h <C 1, both triangular and deformed lattices 
exist only in vicinity of these points. 

For analysis of lattices, it is also useful to derive the 
energy as a function of the average lattice shift, P, and 
the relative phase shift between the neighboring layers, 
4> a . For that we substitute expansion (51 1 up to the linear 



order with the phase gradient given by Eq. ( 52 1 into the 
energy (48 1 and obtain 



£(Va,P) 

L 



-^cos(y> a )sin j cos/3 



h 2 



sin (f a ) 



1 + 2 sin 2 /? cos 2 



hL 



(55) 



In particular, the result ( 53 1 corresponds to the minimum 



of this energy with respect to f a when the condition ( 54 1 



is satisfied. We will see that this relatively simple energy 
function of two variables, whose shape evolves with the 
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FIG. 7: Examples of the energy landscape ( 55 I as a function of 



the average lattice shift (3 and the amplitude of the alternating 
phase deformation (f> a for L = 2 and &/&o — 3.5 and 3. For 
<!> = S^^o the equivalent minima at (<f> a , (3) = (0, 7r) and (tt, 0) 
correspond to the rectangular lattice, while for $ = 3$o the 
minima at (4> a ,P) = (n/2,n/2) and (tt /2,3tt /2) correspond 
to the triangular lattice. 



magnetic field, describes a surprisingly rich behavior in 
the vicinity of the integer-flux-quanta points. The typical 
energy landscapes for the cases of half-integer and integer 
flux quanta per junction are illustrated in Fig. [7] 

Lets study zero-current ground states first. We 
start with stability analysis for the rectangular lat- 
tice which is realized at ip a — 0, tt and gives ground 
state in the most part of phase space. In this case 
the energy is minimal with respect to (3 either at 
(3 = 0, for cos (ip a ) sin(hL/2) > 0, or at — n, for 
cos ((fa) sin (hL/2) < 0, see, e.g., left part of Fig. [7] where 
sin(/iL/2) = — 1. Expanding the energy with respect to 
tp a near the point ip a = 0, 



£(¥>«> o) 5 



S1 " ' 2 J 



hL 



we conclude that the rectangular lattice is stable if 



hL 



> 



L 

7r 



(56) 



which coincides with the condition (54| at /3 = 0,7r. In 



real variables the condition ( 56 1 can be rewritten as 



For the deformed-lattice solution (53 1 the energy at 
fixed (3 is given by 



1 sin 2 if cos 2 (3 
W) ~ Ll+2sin 2 /3cos 2 ^ 



L 



1 + 2 sin 2 /? cos 2 



ML 



(58) 

This energy always has extremums at f3 = 0, n, and tt/2. 
As follows from Eq. (53 1, the state at (3 = tt/2 always 
corresponds to the triangular lattice, (p a = tt/2. We find 
now the conditions that the energy reaches minima at 
these values of (3. Consider first the point (3 = 0. Ex- 
panding the energy near this point, 



„,„s 1 t hL 
£(/3)«- z sin 2 — 



L 



,hL 



L 
h? 

1 + 2 cos 2 



hL 



2L 2 



hL 



W C ° S ~2 



2 \ 2 
we obtain that it corresponds to minimum if 

2i 2 



tan 



hL\ 



1 



W {h 2) 



> 



ti 2 



As L/h <C 1, this inequality is valid almost everywhere 
except in the vicinity of the integer-flux-quanta points 
at which sva(hL/2) — * and cos 2 (hL/2) — > 1 where this 
condition can be rewritten in an approximate simpler 
form, 



hL\ 
1 2 



> 



V2L 
V3h' 



Comparing this condition with the condition ( 56 1 , we can 



see that the deformed lattice gives the energy minimum 
at (3 = only within the narrow region given by 



2 h 
3 < L 



hL\ 



< 1. 



(59) 



In this region the optimum value of <p a for (3 = is given 
by cos{^ a ) = {h/L)\sm{hL /2)\. 

To find if the triangular lattice at the point (3 = tt/2 
gives the local energy minimum, we expand the energy 
(58 1 near this point, (3 = tt/2 — £, 



2tt$ 
$ 



sin 



7T<I> 
$ 



> 



(57) 



Therefore even at small L/h the rectangular lattice is al- 
ways unstable near the integer-flux-quanta points, $ = 
k&o. This is easy to understand: near these points 
the interaction with the boundaries vanishes and even 
small shearing interaction between the neighboring pla- 
nar Josephson-vortex arrays becomes sufficient to induce 
instability with respect to the alternating deformations. 
Further analysis, however, will show that this instability 
takes place when the rectangular lattice does not give al- 
ready the ground state, meaning that the system actually 
experiences a first-order transition. 




, hL 



We can see that the value (3 = tt/2 corresponds to energy 
minimum if 



tan 



2 hL 



< 



2L 2 



l+2cos 2 ^ " h 2 



or, approximately, \s'm(hL/2)\ < \[§Ljh. We can con- 
clude that near the integer-quanta points <f> = fc$o 
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(hL = 2fc7r) the minimum location switches from (3 = 
to (3 = 7r/2. In the intermediate region given approxi- 
mately 



2 /i 
3<L 



~2J 



< V6 



(60) 



the energy has local minimums at both points, (3 = 
and 7r/2. Moreover, in the region (h/L) \sin(hL/2)\ > 1 
the minimum at (3 = is realized by the rectangular 
lattice. This behavior indicates that switching between 
the rectangular and triangular lattices in the ground state 
occurs via a first-order phase transition. 

To find the transition point, we compare the 
triangular-lattice energy, 



/7T 7T 

V2' 2 



L 
h 2 



, hL\ 
2 cos 2 1 



2 ■ 



with the rectangular-lattice energy, 



£ Tect = 8(0,0) = -- 



sin 



hL 



and obtain that the triangular lattice wins if 



hL\ 
1 2 J 



< 



L 

2h. 



1 + 2 cos 2 



hL 
~2 



As this only happens near the points where cos 2 (hL/2) 
1, the equation for the transition points, h t , can again be 
rewritten in a simpler form, 



sin 



htL 



3 L 
2 V 



or in real units, in terms of flux per junction, 

2 



27T$t 







(61) 



(62) 



Comparing Eq. ( |61| with the stability criterion of the 
rectangular lattice ( 56 1 , we indeed can see that before 



the rectangular lattice becomes unstable, it switches to 
the triangular lattice via a first-order phase transition. 
From this equation we can also obtain small shift of 
transition point with respect to the integer-flux-quanta 
point <& = fc<&o as a function of the index k. Writing 
$t = (k + ft.k) $o, we compute 



t,k 



< 1. 



(63) 



The discussed behavior is illustrated in Fig. [8] in which 
the /3-dependences of the energy and current are plotted 
for L = 2 and several values of $ above the point 3<&o- 
The contour plot of energy at the transition point is also 
shown. 
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FIG. 8: The dependences of the energy (left) and current 
(right) on the lattice phase shift (3 for L = 2Aj slightly above 
the point $ = 3$o (curves are marked by the values of 4>/<I>o). 
One can observe several features discussed in details in the 
text. Close the $ = 3$o for $ < 3.1$o the global energy has 
minimum at (3 = tv/2 corresponding to the triangular lattice 
(Eq. (63l actually gives $ t = (3 + l/ir 2 )$o for the transition 



point). Above this point the global minimum is given by the 
rectangular lattice at j3 = 0. Within narrow range of $ the 
energy has minima at both (3 = and 7r/2. The rectangular 
lattice appears at the local minimum at /3 = for $ > 3.08"l?o 
but it becomes unstable with increasing f3. The instability 
points are seen as kinks in the J(f3) curves. In some range of 
$ the absolute value of current has two local maxima within 
< (3 < 7r/2. The critical current switches between these 
maxima with increasing Inset shows the contour plot of 
energy near the transition point, at $ = S.l^o- 



Let us investigate now current-carrying states and be- 
havior of the critical current. In the region given by Eq. 
( 54 1 the current in the rectangular-lattice state is 



J(P) 



hL 



sin/3, 



(64) 



that 



where we assumed for dcfinitcness 
cos ip a sin(/iL/2) > 0. If the condition (54 1 is vio- 



lated then the deformed lattice is realized. In this 
case, differentiating the energy (58 1 with respect to 
(3, we obtain the current in the dt 
(including the triangular lattice) 



eformed-lattice state 



,2 hL 



J(l3): 



(1 + 2 



cos 



2 hL\ 



L (1 + 2 sin 2 cos 2 if Y 



2Lcos 2 

2 sin 2(3, 



h 2 



(65) 



for 



sin 



hL 



cos/3 



Ll + 2sin 2 /?cos 2 ^ 



< 1. 



In particular, at hL = 2irk this formula reproduces re- 
sults ( 46 1 and ( 47 1 obtained from the elliptic- integral rep- 



resentation. 

Consider first the region where the ground state is 
given by the rectangular lattice. Maximum current in 
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this state would be achieved at j3 = ir/2 but the con- 
dition (54 1 always breaks down before that. From this 



condition we compute the value of (3 at which the rect- 
angular lattice becomes unstable 



| cos p t | = - 



2f (l + 2cos 



,2 /iL N 
2 / 



sin 



hL I 



2hL | 8L 2 



_ ro ^( 1+2cos 2M) 
(66) 

In a wide range of parameters, away from the regions 
given by Eq. (61 1, the maximum current is achieved at 



this instability point 

Jc = 



1 2 ) 



sin/3 t 



(67) 



In particular, in most part of the parameter space, for 
| tan(/i£/2)| ^> L/h we obtain much simpler results 



I cos A I 



LI 

h~ 



2 cos 

hL 



2 h_L 
2 



sin 



hL 



L\ 2 (l + 2cos 2 ^)' 
h 



2 sin 



2 hL 
2 



(68) 



(69) 



In this region the critical current is only slightly smaller 
than the "Fraunhofer" result (2/h) |sin (hL/2)\. At 
the half-integer-flux-quanta points hL = (2k + l)ir 
these equations reproduce result ( 42 1 obtained from the 
elliptic-integrals representation. The property that the 
rectangular lattice is always unstable at some lattice dis- 
placement (3 also has important dynamic consequences. 
It means that the lattice can not maintain its static rect- 
angular configuration when it starts to move. 

The critical current has a nontrivial behavior in the 
vicinity of points hL = 2nk where | sin (Ar) | *C 1 and 
general formula ( 65 ) can be simplified as 



2L 



1 



3/i 2 sin 2 (^) \ 

K2J 'sin 2/3. 



h 2 \ 2L 2 (2-cos2f3) , 
This gives the critical current near hL = 27rfc 



J c (h,L) 




2L 2 




for sin(/iL/2) <C L/h. This results shows that the depen- 
dence J c (&) for junction stacks always has local maxima 
at ^ — /c<f>o, in contrast to the Fraunhofer dependence 
for which the critical current vanishes at these points. To 
find the critical current behavior in the whole field range 
in the region h/L 1, we numerically found maximum 
of J(/3) with respect to (3 and different hL = 27r<i>/<i>o 
and L. Figure [9] illustrates the field dependence of crit- 
ical current and current dependence of the lattice struc- 
ture within one oscillation period 2.5<&o < < 3.5<i>o for 
L = 2Xj. To visualize the lattice structures, we represent 
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FIG. 9: The main plot: the field dependence of the critical 
current for L — 2Xj within one oscillation period 2.5$o < $ < 
3.5$o. The grey level codes the value of | cos ip a \ with the light 
grey color in most part of the plot corresponding to the rectan- 
gular lattice | costal = 1 and the black color near O/^o = 3 
corresponding to triangular lattice [cos<^„| = 0. The solid 
line shows the Fraunhofer dependence sin(7r ( I>/$o)/(7r$/<I , o)- 
The inset above shows blowup of the region near $/ ( l > o = 3. 
Plots at the right side illustrate representative dependences 
| J(P)\ for three values of $ marked by arrows. Dashed curves 
in these plots show corresponding dependences for usual small 
junction. The kinks in the J c ($) curve at $/$o ~ 3 ± 0.085 
occur due to the switching between different maxima in the 
|J(/3)| dependence. 



the values of | cos tp a \ by grey level. In most part of the 
current-field diagram the rectangular lattice is realized 
shown by light grey (|cosy> a | = !)■ The triangular lat- 
tice shown by black (| cos ip a \ — 0) appears in the ground 
state only in vicinity of the point <!> = 3<&o. Exactly at 
this point the lattice remains triangular up to the crit- 
ical current. Slightly away from this point the lattice 
deforms with increasing current. In the range of param- 
eters given by Eq. (60 1 the dependence |J(/?)| has two 
maxima within < j3 < tt/2 (see left plot in Fig. [8]). As 



a consequence, the field dependence of the critical current 
has kinks related to switching between these maxima. 



VII. SLOW DYNAMICS IN OVERDAMPED 
REGIME: OSCILLATIONS OF THE FLUX-FLOW 
VOLTAGE 

When the external current flowing across the layers 
exceeds the critical current, the lattice starts to move. 
In general, dynamic behavior is quite complicated. A 
simple situation is realized only at slow lattice motion 
in the overdamped case when the lattice deformations 
have time to adjust to the current lattice position. In 
this case the lattice moves in the periodic potential given 
by its static energy (12 1 and one can use static results 



to predict the I-V dependences. On the other hand, one 
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FIG. 10: The representative field dependences of mean- 
squared average of current, (J 2 (a)}, with respect to the lat- 
tice displacements which determines the amplitude of relative 



voltage oscillations SU/Uff at slow velocities via Eq. (72 l 
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FIG. 11: The field dependence of the ratio of voltage- 
oscillation maxima at integer-flux-quanta points (5U max> 2) 
and at half-integer-flux-quanta points (5U m ax,i)- The in- 
set illustrates definitions of dU m ax,i and 8U m ax,2 in the 
schematic voltage-field dependence. For comparison we also 
show plot of the ratio of the critical-current maxima squared 
( Jmax,2/ Jmax,i) 2 ■ Extraction of the typical field B 1 / 2 from 
the analysis of the voltage oscillations allows accurately eval- 
uate the anisotropy factor 7 using Eq. ( 73 1 . 



can expect that the voltage oscillations are less sensitive 
to inhomogcncitics than the critical-current oscillations, 
because the homogeneously moving lattice smears away 
disorder. The critical current is also smeared by thermal 
fluctuations. These are possible reasons why it is easier 
to observe and interpret the magnetic oscillati ons in the 
flux-flow resistivity than in the critical current. ^ * 15 * 16 * 17 ! 

In general, the dynamic behavior also depends on the 
dissipation mechanism. In BSCCO in a wide range of 
magnetic fields the flux-flow resistivity is mainly deter- 



mined by the in-plane quasiparticle conductivity a a i 
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Only wh en the magnetic field exceeds a typical value 
B a = \J Oab I &c®o I (V27T7 2 s 2 ) , the c-axis conductivity, 
<7 C , gives dominating contribution to the flux-flow dissipa- 
tion. In this limit the flux-flow resistivity becomes field- 
independent. In BSCCO the field B a is typically several 
times larger than the crossover field B cr . An important 
feature of the in-plane dissipation regime at B < B a is 
that the lattice velocity at fixed applied current is very 
sensitive to lattice structure, the smallest velocity is re- 
alized for the triangular lattice and the largest velocity is 
realized for the rectangular lattice. As the lattice struc- 
ture in the regime B > Bl continuously changes with lat- 
tice displacement, the dynamic behavior in this regime is 
rather complicated. To avoid this complications, we limit 
ourself here by a simple case of dominating c-axis dissi- 
pation in the crossover region, Bl > B a . In this case 
the viscous-friction coefficient weakly depends on lattice 
structure. 

In the case of structure-independent viscous-friction 
coefficient Vff, time variation of the lattice phase shift 
obeys equation 



d& ~j / \ 
1/ ff- (it + J ( a ) = J c 



(70) 



where J cx t is the external current, the current J{a) = 
J(a, hL, h) is given by Eq. ( 13 1 (for brevity we again 
skip in equations its dependence on the magnetic field 
and size), and the viscosity coefficient, Vff, is related to 
the flux-flow resistance of the stack, Rff, 



v ff 



N<S>o 



2ncR 



ff 



where N is the number of junctions in the stack. The 
voltage drop per one junction U is related to da/dt by 
the Josephson relation 



U = 



<I>o da 

2-7TC dt 



Solution of Eq. ( 70 ) is given by the implicit relation 

f ^f da' 
Jo Jext ~ J (a 1 ) 

from which we obtain the average phase change rate, 

n -1 



da 
~dt 



vffda 



1 J ext 

and the flux-flow voltage 



U 

Wf 



j, 



ext 
7T 



J(a) 



da 



Jext - J (a) 



(71) 



with Uff = Rff J being the bare flux-flow voltage with- 
out the periodic potential. As the current J(a) = 
J(a, hL, h) oscillates with the magnetic field, this flux- 
flow voltage will also experience similar field oscilla- 
tions. In particular, when the external current signifi- 
cantly exceeds the critical current, J cxt S> J c , we ob- 
tain weak oscillating correction to the flux-flow voltage, 
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8U = U-U„, 



8U/U ff n-{j\a))/Jl 



(72) 



where (/(a)) = (1/tt) f£ f(a)da is the average with re- 
spect to the lattice phase shift. As SU oc ^J 2 (a)), the 
behavior of SU is overall similar to the behavior of the 
critical current but the amplitude of voltage oscillations 
roughly scales as the critical current squared. Figure [T0| 
shows the field dependences of the average (J 2 (a)}, which 
determines the amplitude of weak voltage oscillations, for 
three junctions sizes, L = 3, 4, and 6. 

Consider in more details behavior of SU at the points 
<f> = (j/2)<I>o. To find the amplitude of the small volt- 
age correction, SU ma x.i, at the half-integer flux quanta 
points, <f> = (k+ 1 /2)<E»o , we have to find the a-average of 
Jf(ct) where the current J\(a) is given by Eq. (41 1 with 
the parameter rh given by Eq. ( 39 1 . Similarly, the am- 



plitude of the voltage oscillation at the points $ = k$o, 
which we notate as SU max _2, is determined by (JfC 01 )) 
where the current J2(«) is given by Eq. (45 1 with the 
parameter m given by Eq. ( 43 ) . 



Figure 
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shows the 

computed field dependence of the ratio SU max ^.l SU max .\ 
together with the ratio ( Jmax.2/ Jmax,i) 2 ■ In practice, 
to extract the ratio SU max _2/ SU max ,i from experimental 
voltage-field dependence, one should plot smooth curves 
via local maxima and two sets of local minima as it is il- 
lustrated in the inset of Fig. |11| subtract the two minima 
curves from the maxima curve, and compute the ratio of 
the differences. One can see that the field dependences of 
SU max ,2/SU m axA and ( J ma x,2/ J m ax,i) 2 are almost identi- 
cal. This plot gives possibility for accurate determination 
of the anisotropy factor from the voltage oscillations us- 
ing the field scale. In particular, the ratio SU2/SU1 drops 
to 0.5 at the field B 1/2 = 1.302B L = 1.302<S> L/(2nj 2 s 3 ) 
(see Fig. 11 I meaning that 7 can be extracted from this 
field as 



7 w 330. 



' L[fj,m] 

£1/2 py 



(73) 



where we used the BSCCO interlayer spacing s w 1.56 
nm. For example, using data reported in Ref. [15] for the 
sample H55 with L = 5.5 /im, we estimate -B1/2 ~ 3.6T 
giving a very reasonable estimate 7 s» 408. This estimate 
is significantly larger than the value 7 ~ 110 obtained by 
the authors themselves and the difference comes from the 
value of the numerical constant in the crossover field. 

We have to mention that an alternative mechanism of 
<i>o/2-periodic voltage oscillations at fixed current exists 
at high lattice velocities due to switching between the 
Fiske steps. 26 However, experimentally, the most regular 
voltage oscillations are obse rved at vo ltages much smaller 
than the first Fiske volt age piaiHE] 



VIII. SUMMARY 

In summary, we studied magnetic oscillations of the 
critical current and lattice configurations in stacks of in- 



trinsic Josephson junctions, which are realized in mesas 
fabricated from layered high-temperature superconduc- 
tors. Depending on the stack lateral size, oscillations may 
have either the period of half flux quantum per junction 
(wide-stack regime) or one flux quantum per junction 
(narrow-stack regime). We studied in detail the crossover 
between these two regimes. Typical size separating the 
regimes is proportional to the magnetic field meaning 
that the crossover can be driven by the magnetic field. In 
the narrow-stack regime the lattice structure experiences 
periodic series of first-order phase transitions between 
aligned rectangular configuration and triangular config- 
uration. The triangular configurations in this regime is 
realized only in narrow regions near magnetic-field val- 
ues corresponding to integer number of flux quanta per 
junction. For slow lattice motion similar crossover can 
also be observed in the oscillations of the flux-flow resis- 
tivity. Quantitative study of the crossover allows for a 
very accurate evaluation of the anisotropy factor. 
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APPENDIX A: REGIONS OF MONOTONIC 
SAME-SIGN SOLUTIONS FOR v(u) 

One can distinguish two types of monotonic solu- 
tions depending on wether or not the smooth phase 
v(u) = (tp(u) — 7r/2)/2 changes sign inside the junction 
(see Fig. |2j. For the changing-sign solution the condi- 
tion < m < 1 always holds. In this Appendix we find the 
boundary values of a separating these two types of mono- 
tonic solutions. For definiteness, we consider the region 
2kn<hL<(2k+l)n and -n/2+hL-2kir<a<ir/2 (grey 
area in the lower part of the phase diagram in Fig. |2j. 
For the monotonic changing-sign solution in this range 
we have 0<ipo<ir/2, 7r/2<0£<7r, i.e., 



ipo = arcsm 



<Pl = 7r — arcsm 




l + 2cos 2 (/iL~a) 



There are two boundaries in the region, one correspond- 
ing to the condition 0o = tt/2 (uq = 0) below a = tt/2 
and another corresponding to 4>l = ir/2 (vl = 0) above 
a = —7T /2+hL~2kTT (see Fig. [2]). For the first boundary, 
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ao(h,L), from Eqs. (20al and (22) we obtain 



/mo 



m 



l/m 



K(mn)— F ( arcsin 4 / — 

1 07 I Y l + 2cos 2 (/iL-a ) 



1 



rn 



1 + 2 cos 2 (a ) 



(Ala) 
(Alb) 



where if (to) = F(tt/2, m) is the complete elliptic in- 
tegral of the first kincH^J and we used the identity 
F (tt — j3, m)—K(m) — K(m)—F ((3, m). Analyzing similar 
equation for the second boundary aj,(/i, L), we find that 
it is related to ao(^) L) as ai(/i, L) = hL—2kTr—ao(h, L). 
One can check that a>o(h, L) — > 7r/2 for ftX — > 2fc7r, (2fc + 
l)vr and for L — > oo, i.e., the region of the same-sign 
monotonic solution vanish in these limits. 



One can neglect shift of rj due to the finite-size correc- 
tions to vq and vl- Without interaction between edge 
deformations, the energy can be written asP 

£o = — sin (2t>oo) cos (a) — — sin {2vlo) cos (hL — a) 
h h 

+ U2° dvVl ~ cos{4v) ~ L ^ 

= — ^ (2-v / 2+cos2a-j2+cos2(a-hL))-^r. 
hy/2 V v V h? 

where «oo and vlq are the surface deformations neglecting 
finite-size correction, 



tan vqo 
tan vlo 



l-VT+2 



cos z a 



•\/2cosa 
1 - v/l + 2cos 2 (a-hL) 
•\/2cos (a — hL) 



APPENDIX B: WEAK FINITE-SIZE EFFECTS 
AT LARGE L 

In this Appendix we derive finite-size corrections to 
the critical current due to interaction between the surface 
solitons. Consider for definiteness the case of monotonic 
solution. The nonmonotonic solution can be treated sim- 



ilarly. In the limit L/h^> 1 the parameter to in Eq. (21 ) 
is close to one. Separating small correction, to = 1 — r}/2 
with rj <c 1, we evaluate the integral as 



dv 1/32 

— = ss - m tan vn tan vl 

Vl + ^/2-cos 2 (2v) 2 V V 



This gives the following result for rj 



-32 tan vq tan vl exp (—"v8L/h 



(Bl) 



corresponding to the elliptic-function parameter to S3 
1 + 16 tan vq tan^i exp (— \/%L/ti\ . The boundary con- 
ditions can be represented as 



cos (2vq) 
cos (2vl) 



l + 7?/4 



v/l + 2cos 2 (a)' 

l+ji/4 
v /l+2cos 2 (/ii-a)' 



The finite-size correction to the energy change can now 
be estimated as 



1 /-vlo , > 

- / dv ^v/l + ry-cos (4^)-^! -cos (4?;) J 



t'oo 

V 



S£ = 

~ /i%/32 

Therefore, the total total energy can be written as 
L 



Li] 
2h? 



£(a, h, hL) 



h 2 



+ -j=- (2 - y/2 + cos(2a) - V2T cos [2 (a — hL)] 

8\/2 cos(a) cos (a — hL) exp (—y/8L/ti) 
h (1 + ^2+ cos (2a)) (l + ^2+cos[2(a-/iL)]) . 

(B2) 



The last correction term describes the exponentially 
small interaction energy between the surface solitons. It 
is important to note that the finite-size correction breaks 
7r-periodicity with respect to parameter hL meaning that 
the states with $ = fc<& and $ = (k + l/2)$ are not 
equivalent any more. 

For the Josephson current we obtain 



J(a, h, hL) = 



sin 2a 



sin 2 (a — hL) 



V2h \ ^2 + 008 20 y/2 + cos 2 (a- hL) 



8V2 



exp 



VBL\ d_ 
h I da 



cos a cos (a — hL) 



1 + V2 + cos 2a) (l + ^J2 + cos 2 (a - hL) 



Assuming that without the finite-size correction the maximum current flows at a = a m {hL), we obtain for the 
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finite-size correction to the critical current 



oJ c (h,hL) = — - — exp 



h 



h I da 



cos a cos (a — /iL) 



(l + y/2 + cos(2a)) (1 + ^2 + cos [2 (a — hL)] 



OL—ot m {hL) 



In particular, near the maxima hL = irj (<f> = j<f>o/2), using the resuliPa m (0) = 0.921, we obtain result ((37}. 
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